/*
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef FRACTALNOISEGENERATOR_H_
#define FRACTALNOISEGENERATOR_H_

#include "OgreMath/OgreMath.h"
#include <stdlib.h>

class FractalNoiseGenerator
{
private:
	static float***
	allocMIPMaps(unsigned int iResolution)
	{
		float*** pMaps=new float**[iResolution];

		for(unsigned int iMip=0; iMip<iResolution; iMip++)
		{
			const unsigned int IMGSIZE=(1<<(iResolution-iMip-1));
			pMaps[iMip]=new float*[IMGSIZE];

			for(unsigned int x=0; x<IMGSIZE; x++)
				pMaps[iMip][x]=new float[IMGSIZE];
		}
		return pMaps;
	}

	static void
	freeMIPMaps(unsigned int iResolution, float*** pMaps)
	{
		for(unsigned int iMip=0; iMip<iResolution; iMip++)
		{
			const unsigned int IMGSIZE=(1<<(iResolution-iMip-1));
			for(unsigned int x=0; x<IMGSIZE; x++)
				delete pMaps[iMip][x];

			delete pMaps[iMip];
		}
		delete pMaps;
	}

	static float**
	allocHeightMap(unsigned int iResolution)
	{
		const int IMGSIZE=(1<<iResolution);
		float** pMap=new float*[IMGSIZE+1];

		for(int x=0; x<IMGSIZE+1; x++)
			pMap[x]=new float[IMGSIZE+1];
		return pMap;
	}

	static void
	freeHeightMap(unsigned int iResolution, float** pMap)
	{
		const int IMGSIZE=(1<<iResolution);

		for(int x=0; x<IMGSIZE+1; x++)
			delete pMap[x];
		delete pMap;
	}
public:
	/*
	* generate 2D fractal noise using square diamond algorithm
	* @return float[(1<<iResolution)+1][(1<<iResolution)+1] array with noise values
	*/
	static void
	generate(int iResolution, float fRoughness, float** ppfFNA)
	{
		int IMGSIZE=(1<<iResolution);

		ppfFNA[0][0]=0;
		ppfFNA[IMGSIZE][0]=0;
		ppfFNA[0][IMGSIZE]=0;
		ppfFNA[IMGSIZE][IMGSIZE]=0;

		int iStep=IMGSIZE;
		float fCurrRoughness=1;
		while(iStep>1)
		{
			fCurrRoughness*=fRoughness;
			for(int x=0; x<IMGSIZE; x+=iStep)
			{
				for(int y=0; y<IMGSIZE; y+=iStep)
				{
					//square step
					float fOffset=Ogre::Math::RangeRandom(-1.0f, 1.0f)*0.5f*fCurrRoughness*(float)iStep/(float)IMGSIZE;
					//calc average of surrounding square
					float fSum=0;
					fSum+=ppfFNA[x][y];
					fSum+=ppfFNA[x+iStep][y];
					fSum+=ppfFNA[x][y+iStep];
					fSum+=ppfFNA[x+iStep][y+iStep];
					//set center point of square
					ppfFNA[x+(iStep>>1)][y+(iStep>>1)]=fOffset+fSum/4.0f;
					//a little hack right here: since this generator is used to create a cutting plane
					//and this very point is the origin of that plane, we thereby make sure,
					//that the cutting plane will not exceed the original geometry
					if(iStep==IMGSIZE)
						ppfFNA[x+(iStep>>1)][y+(iStep>>1)]=0;

					//diamond step
					ppfFNA[x+(iStep>>1)][y]=(ppfFNA[x][y]+ppfFNA[x+iStep][y])/2.0f+
							Ogre::Math::RangeRandom(-1.0f, 1.0f)*0.5f*fCurrRoughness*(float)iStep/(float)IMGSIZE;
					ppfFNA[x][y+(iStep>>1)]=(ppfFNA[x][y]+ppfFNA[x][y+iStep])/2.0f+
							Ogre::Math::RangeRandom(-1.0f, 1.0f)*0.5f*fCurrRoughness*(float)iStep/(float)IMGSIZE;
					ppfFNA[x+(iStep>>1)][y+iStep]=(ppfFNA[x][y+iStep]+ppfFNA[x+iStep][y+iStep])/2.0f+
							Ogre::Math::RangeRandom(-1.0f, 1.0f)*0.5f*fCurrRoughness*(float)iStep/(float)IMGSIZE;
					ppfFNA[x+iStep][y+(iStep>>1)]=(ppfFNA[x+iStep][y]+ppfFNA[x+iStep][y+iStep])/2.0f+
							Ogre::Math::RangeRandom(-1.0f, 1.0f)*0.5f*fCurrRoughness*(float)iStep/(float)IMGSIZE;
				}
			}
			iStep>>=1;
		}
	}

	static void
	generateMIPs(int iResolution, float** ppfFNA, float*** pMinMIPs, float*** pMaxMIPs)
	{
		float** pMinSource=nullptr;
		float** pMaxSource=nullptr;

		for(int iMip=0; iMip<iResolution; iMip++)
		{
			if(!iMip)
			{
				pMinSource=ppfFNA;
				pMaxSource=ppfFNA;
			}
			else
			{
				pMinSource=pMinMIPs[iMip-1];
				pMaxSource=pMaxMIPs[iMip-1];
			}

			const unsigned int IMGSIZE=(1<<(iResolution-iMip));

			float** pMinTarget=pMinMIPs[iMip];
			float** pMaxTarget=pMaxMIPs[iMip];

			for(unsigned int iMIPx=0; iMIPx < (IMGSIZE>>1); iMIPx++)
			{
				float* pMinTargetRow=pMinTarget[iMIPx];
				float* pMaxTargetRow=pMaxTarget[iMIPx];

				for(unsigned int iMIPy=0; iMIPy < (IMGSIZE>>1); iMIPy++)
				{
					pMinTargetRow[iMIPy]=Ogre::Math::POS_INFINITY;
					pMaxTargetRow[iMIPy]=Ogre::Math::NEG_INFINITY;
				}
			}

			for(unsigned int x=0; x<IMGSIZE; x++)
			{
				int iMIPx=x>>1;

				float* pMinTargetRow=pMinTarget[iMIPx];
				float* pMaxTargetRow=pMaxTarget[iMIPx];

				float* pMinSourceRow=pMinSource[x];
				float* pMaxSourceRow=pMaxSource[x];

				for(unsigned int y=0; y<IMGSIZE; y++)
				{
					int iMIPy=y>>1;

					if(pMinSourceRow[y] < pMinTargetRow[iMIPy])
						pMinTargetRow[iMIPy] = pMinSourceRow[y];

					if(pMaxSourceRow[y] > pMaxTargetRow[iMIPy])
						pMaxTargetRow[iMIPy] = pMaxSourceRow[y];
				}
			}
		}
	}

	class Constructor
	{
	private:
		float **heightMap;
		//[mip][x][y]
		float ***maxMIPs;
		float ***minMIPs;

		const int expo;
		const int zRange;
		const bool flat;

		const int imgSize;
		const int halfSize;

	public:
		Constructor(const Constructor &rhs) :
			expo(rhs.expo), zRange(rhs.zRange), flat(rhs.flat),
			imgSize(rhs.imgSize), halfSize(rhs.halfSize)
		{
			heightMap = allocHeightMap(expo);
			for (int x = 0; x < imgSize+1; x++)
				memcpy(heightMap[x], rhs.heightMap[x], (imgSize+1)*sizeof(float));

			maxMIPs=allocMIPMaps(expo);
			minMIPs=allocMIPMaps(expo);

			generateMIPs(expo, heightMap, minMIPs, maxMIPs);
		}
		/*inline Constructor& operator = ( const Constructor& rhs )
		{
			expo = rhs.expo;
			zRange = rhs.zRange;
			imgSize = rhs.imgSize;
			halfSize = rhs.halfSize;
			flat = rhs.flat;

			heightMap = allocHeightMap(expo);
			for (int x = 0; x < imgSize+1; x++)
				memcpy(heightMap[x], rhs.heightMap[x], (imgSize+1)*sizeof(float));

			maxMIPs=allocMIPMaps(expo);
			minMIPs=allocMIPMaps(expo);

			generateMIPs(expo, heightMap, minMIPs, maxMIPs);

			return *this;
		}*/
		Constructor(int _expo, int _zRange, float roughness, bool _flat) :
			expo(_expo), zRange(_zRange), flat(_flat),
			imgSize(1<<expo), halfSize(imgSize>>1)
		{
			heightMap = allocHeightMap(expo);
			generate(expo, roughness, heightMap);

			float noiseMax = 0.0f;
			for (int x = 0; x < imgSize; x++)
				for (int y = 0; y < imgSize; y++)
					if (heightMap[x][y] > noiseMax) noiseMax = heightMap[x][y];

			float multiplier = (float)zRange/(noiseMax+0.001f);
			for (int x = 0; x < imgSize; x++)
				for (int y = 0; y < imgSize; y++)
					heightMap[x][y] *= multiplier;

			//generate MIPs
			maxMIPs=allocMIPMaps(expo);
			minMIPs=allocMIPMaps(expo);

			generateMIPs(expo, heightMap, minMIPs, maxMIPs);
		}
		~Constructor()
		{
			freeHeightMap(expo, heightMap);
			freeMIPMaps(expo, maxMIPs);
			freeMIPMaps(expo, minMIPs);
		}

		const int getSize() const {return expo;}

		const VoxelType fn(const Box<Vector3i>& box) const
		{
			const int bottom = (-imgSize)-zRange;
			const int minX = box.min.x + halfSize;
			const int minY = box.min.y + halfSize;
			const int maxX = box.max.x + halfSize;
			const int maxY = box.max.y + halfSize;

			if(!flat)
			{
				//exclude outside areas
				if (maxX < 0 || minX >= imgSize ||
					maxY < 0 || minY >= imgSize ||
					box.min.z > zRange || box.max.z < bottom)
					return Nil;

				if (box.max == box.min)
				{
					if (box.min.z <= heightMap[minX][minY])
						return Solid;
					return Nil;
				}

				//create block beneath the noise
				if(box.max.z <= -zRange &&
					box.max.z > bottom)
				{
					if((box.max-box.min).maxComponent()<10)//not too coarse
						return Solid;
					else
						return Both;
				}

				if(minX < 0 || maxX >= imgSize ||
						minY < 0 || maxY >= imgSize)
					return Both;
				//convert box to mip pixel
				Box<Vector2i> rect( Vector2i(minX, minY),
						Vector2i(maxX, maxY) );

				int iMIPExpo=AuxMath::log2(std::max(maxX-minX, maxY-minY) );

				int iCoordExpo=iMIPExpo+1;

				Vector2i vMin=Vector2i(minX>>iCoordExpo, minY>>iCoordExpo);

				Vector2i vMax=Vector2i(maxX>>iCoordExpo,
										maxY>>iCoordExpo);
				//fail-safe
				if(iMIPExpo>=expo)
				{
					iMIPExpo=expo-1;
					vMin=vMax=Vector2i(0);
				}

				float fMin=minMIPs[iMIPExpo][vMin.x][vMin.y];
				float fMax=maxMIPs[iMIPExpo][vMin.x][vMin.y];

				if(vMin!=vMax)
				{
					Vector2i vDiff=vMax-vMin;
					//sample pixels
					for(unsigned int x=0; (int) x <= vDiff.x; x++)
						for(unsigned int y=0; (int) y <= vDiff.y; y++)
						{
							if( (!x && !y))
								continue;

							float fMinSample=minMIPs[iMIPExpo][vMin.x+x][vMin.y+y];
							float fMaxSample=maxMIPs[iMIPExpo][vMin.x+x][vMin.y+y];

							if(fMinSample<fMin)
								fMin=fMinSample;
							if(fMaxSample>fMax)
								fMax=fMaxSample;
						}
				}

				if (box.max.z <= fMin)
					return Solid;
				if (box.min.z > fMax)
					return Nil;
			}
			else
			{
				if (maxX < 0 || minX >= imgSize ||
					maxY < 0 || minY >= imgSize ||
					box.min.z > zRange)
					return Nil;

				if (box.max == box.min)
				{
					if (box.min.z <= heightMap[minX][minY])
					{
						float minNeighbor = heightMap[minX][minY];
						if (minX > 0)
						{
							minNeighbor = std::min(minNeighbor, heightMap[minX-1][minY]);
							if (minY > 0) minNeighbor = std::min(minNeighbor, heightMap[minX-1][minY-1]);
							if (minY < imgSize-1) minNeighbor = std::min(minNeighbor, heightMap[minX-1][minY+1]);
						}
						if (minX < imgSize-1)
						{
							minNeighbor = std::min(minNeighbor, heightMap[minX+1][minY]);
							if (minY > 0) minNeighbor = std::min(minNeighbor, heightMap[minX+1][minY-1]);
							if (minY < imgSize) minNeighbor = std::min(minNeighbor, heightMap[minX+1][minY+1]);
						}
						if (minY > 0) minNeighbor = std::min(minNeighbor, heightMap[minX][minY-1]);
						if (minY < imgSize-1) minNeighbor = std::min(minNeighbor, heightMap[minX][minY+1]);
						if (box.min.z >= static_cast<int>(minNeighbor)) return Solid;
					}
					return Nil;
				}
				return Both;
			}
			return Both;
		}
	};
};

#endif /* FRACTALNOISEGENERATOR_H_ */
